function vOut=compute_product_QSE(G,vIn)
%compute_product_QSE -- matrix-vector product in QSEPACK.
%   Given a metric structure G returned by initialize_metric_QSE(...) 
%   and a column vector vIn containing algebraic state coordinates, then 
%
%      vOut = compute_product_QSE(G,vIn);
%
%   returns their matrix-vector product in the column vector vOut.
%
%   The indexing convention for flattening the three xi-coefficient 
%   indices {j,k,l}
%
%      xi{j,k}(1,l) = ${}\sp{j}\xi^{k}_{l}$ (in LaTeX notation)
%
%   into a single index for vIn and vOut is as follows:
%
%     index j varies most slowly, over [1:G.rank].
%     index k varies next most slowly, over [1:G.order], and
%     index l varies most rapidly, over [1:G.lengths(k)],
%
%   or in MATLAB notation (assuming xi is a cell array of row vectors)
%
%   vIn = reshape(cell2mat(xi).',[],1)

% ******************************************************************************
% Copyright (C) 2008 John Arthur Sidles.  This copyright is assigned  
% to the Institute for Soldier Healing of Seattle, Washington, USA.
% ------------------------------------------------------------------------------
% This source file is part of the QSEPACK Template Library (release 1.0d).  
% The QSEPACK Template Library is free software: you can redistribute it 
% and/or modify it under the terms of the GNU General Public License 
% as published by the Free Software Foundation, either version 3 of 
% the License, or (at your option) any later version. 
% 
% The QSEPACK Template Library is distributed in the hope that it will  
% be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
% GNU General Public License for more details. 
% 
% You should have received a copy of the GNU General Public License along  
% with the QSE Template Library. If not, see <http://www.gnu.org/licenses/>.
% 
% The author, John A. Sidles, can be contacted through the Institute for 
% Soldier Healing, 4532 48th Av NE, Seattle, Washington, 98105 USA, or
% by email at "jsidles@soldierhealing.org"
% ------------------------------------------------------------------------------
% The QSEPACK template library is maintained using the open-source
% literate programming tool nuweb (v1.1.1), by Preston Briggs.
% ------------------------------------------------------------------------------
% The QSEPACK template library is distributed as a gzipped tarball file 
% "QSEPACK_v1_0d.tar.gz" that includes the master documentation file 
% "QSEPACK_v1_0d.pdf". 
% ******************************************************************************

% ------------------------------------------------------------------------------
% scrap: code for compute_product_QSE.m
%   G.rank    = [1] number of block rows and columns (of $G$-blocks in Table 2D)
%   G.cvecs   = {rank,rank}[ncoef,1] cell array of vectors ($c$ in Table 2E)
%   G.smats   = {rank,rank}[ncoef,ncoef] cell array of sparse matrices ($B+D$ in Table 2F)
vIn  = reshape(vIn,[],G.rank);     % reshape column vector into blocks
vOut = complex(zeros(size(vIn)));  % preallocate output storage
for i1 = 1:G.rank      % loop over metric block rows
  for i2 = 1:G.rank    % loop over metric block columns
    vOut(:,i1) = vOut(:,i1) ...    % efficiently compute-and-increment product
      + G.cvecs{i1,i2}*((G.cvecs{i2,i1}')*vIn(:,i2)) ...
      + G.smats{i1,i2}*vIn(:,i2);
  end
end
vOut = reshape(vOut,[],1); % restore to column vector output